% Main program for dynamic phase screen generation using Zernike method
% Written by Benchu lu; email:benchuul@163.com; Taiyuan University of Technology
% Reference: Benchu Lu and Dongmei Cai, "Simulation of a temporally dynamic
% atmospheric wavefront using the modified von-Kármán spectrum and 
% Zernike-based method," Opt. Express 32, 38609-38629 (2024)
clear; close all; clc;
% Initializing the parameters for dynamic phase screen
%------------------------------------------------------------
params.n_max =25; % maximum n value for Znk Plyns
params.j_max = Noll_nm_to_j(params.n_max, params.n_max)+1;
params.D = 2; % aperture diameter [m]
params.resol = 256; % resolution of phase screen [pixels per side]
params.lambda = 525e-9; % wavelength [m]
params.L0 = inf; % Outer sacle [m]
params.l0 = 0; % Inner scale [m]
% sp:1 for spherical wave; 0 for plane wave；当sp=1时，计算球面波要求multilayer=1
params.sp = 0; 
params.Cn2 = 1e-16; % refractive-index structure parameter [m^(-2/3)]
params.L = 20e3; % propagation distance [m]
% Fried parameter r0 [m]: r0可由Cn2(z),L和lambda得出
params.r0 = r0_generate(params.Cn2, params.lambda, 0, params.L, params.sp, params.L); 
params.W_V = 10; % Wind Velocity [m/s]
params.dynatime = 2; % dynamic time for phase screen [s]
% 时间序列长度 dynatime*samplehz 总为向下取整的偶数，有单数时间采样频率需求的可自行修改代码
params.samplehz = 500;  %sampling frequence for 1 s [Hz],
% RotateAgl: 风速方向相对于x轴的夹角   0：x方向； pi/2：y方向
params.RotateAgl = (1/4) * pi; % [rad]
%----------------------------------------------------------------
%------------------分层廓线模型参数设置---------------------------
params.multilayer = 0; % 1 for multi-layer case; 0 for homogeneous medium
params.layerNum = 5; % 多层大气模型的等间距分层数
% W = 21, A = 1.7e-14, r0 = 0.05, theta_0 = 6.9e-6
% W = 10.5, A = 0.25e-14, r0 = 0.10, theta_0 = 12.9e-6
%params.W = 10.5; params.A = 0.25e-14;
%params.Cn2 = HVB_Cn2(params);
%params.W_V = Bufton_W_V(); % Bufton风速模型 H_V(z) :z~[km]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算系数时间功率谱，协方差矩阵和K-L函数
[W_aj, C, K_L, x, idx, ph, params] = ActionFun(params);
disp('The temporal spectrum and K-L functions have been computed');
% end of the main program